
      !<@file
      !>@brief Rotina para geração dos sismogramas.
      !>
      !> permite a geração dos sismogramas através da solução numérica da EDP (no domínio do tempo ou da frequência), por meio de métodos numéricos (MDF) .
      !>@todo implementar fontes e receptores fora do nó da malha
      !>@todo implementar outro esquema para relacionar fontes e receptores para facilitar a entrada dos casos "reais" 1 fonte e receptores em todo o cabo.... outra fonte, receptores em todo cabo
      !>@todo lembrar dos parametros de aquisicao: amostragem 2milisegundos, tempo de gravaçao 9segundos, distancia da fonte para o primeiro receptor 45m, distancia entre receptores 12.5m, tamanho do cabo de 3 a 9km, o navio se desloca 25m para o proximo tiro
      subroutine simulacao(tsim,filesism,iprint)

      use globais
      implicit none

!	inversao 
!       double precision, dimension(:), allocatable :: X !< variaveis de otimizacao/parametros da inversao

!       double precision:: cmin0,cmax0
!       integer :: m !< numero de medidas 
      integer :: i,j
!       double precision, dimension(:), allocatable :: Xaux !< variaveis de otimizacao/parametros da inversao auxiliar.


      double precision:: dreceiv !< distância entre receptores, variável auxiliar na geração dos sismogramas
!       double precision, dimension(:), allocatable :: Gaux
!       double precision :: Faux

!       double precision :: VALUE	
!       integer :: kncl
!       integer :: knclbest !< Melhor número de classes para clusterização.
!       doubleprecision :: valuebest !< Melhor valor da função objetivo entre os números de classe testados


! 
! 	integer :: MaxIter
! 
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! ! non linear conjugate gradient, napack, netlib
! ! 	!input
! 	double precision :: STEP  !--STEP SIZE ALONG SEARCH DIRECTION FOR FINAL ITERATION		   
!  	double precision :: Tol  !--COMPUTING TOLERANCE (ITERATIONS STOP WHEN MAX-NORM OF GRADIENT .LE. Tol)  
! 	integer :: limit
! 	double precision , dimension(:,:), allocatable::Work
! 	!output
! 	double precision :: E  !--MAX-NORM OF GRADIENT

! 
! 	EXTERNAL BOTH,GRAD,PRE,VALUE
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! 
! 
!       !< L-BFGS 
! !       PARAMETER(NDIM=2000,MSAVE=7,NWORK=NDIM*(2*MSAVE +1)+2*MSAVE)
!       DOUBLE PRECISION, dimension(:), allocatable:: DIAG
!       DOUBLE PRECISION, dimension(:), allocatable:: W
!       DOUBLE PRECISION :: EPS,XTOL,GTOL,STPMIN,STPMAX
!       INTEGER ::  IPRINT(2),IFLAGlbfgs,ICALL,MP,LP
!       LOGICAL ::  DIAGCO
!       integer :: msave
!       integer :: icallmax

! 
!       EXTERNAL LB2
!       COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
! 
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 



      integer :: tsim!< tipo de simulação tsim=1 MDF, domínio da frequência; tsim=2 MDF, domínio do tempo


      ! para simulação no domínio da frequencia
      double precision, dimension(:,:), allocatable :: P1,P2,P3
      double precision, dimension(:,:), allocatable :: Caux,Aaux
      double precision :: cminmodel,cmaxmodel ! velocidades mínima e máxima no modelo
      double precision :: deltatmax !< valor máximo para delta t
      double precision :: deltamax !< valor máximo para delta 
      double precision :: deltat !< passo de tempo

      integer :: plotsnap !< plotsnap = 1 salva o campo em todo domínio em bin para futuro snapshot
    

      integer :: ierr

      double precision :: lbeg,lend,ldnum

      double precision :: t ! provisorio     
       integer :: it ! provisorio 
       double precision :: dt ! para sismograma e snapshots     

!       integer :: simulacao
      character*10 :: filename4
      character*5 :: filename1,filename2,filename3
      character*1 :: barra= achar(92),pl= achar(39),plpl= achar(34)     
      character*2 :: charit

!       double precision :: mapCemX

      complex*16 ::Ptrecaux
!       double precision :: residuo
      character*5 :: chks
      character*2 :: chigsr

      real, dimension(:,:), allocatable:: PtrecRaux !< Copia dos sismogramas.
      integer :: it1,it2,itIN,nttotalIN,ntpsIN
      double precision :: t1,t2,tIN,deltaTIN

      character*5 ::filesism !< primeiros 5 caracteres do nome dos arquivos referentes ao sismograma
      integer :: iprint !< indica se plotará na tela as chamadas ao pswigp

      write(*,*) '# simulacao'




      write(*,*) '#menor comprimento de onda (no modelo) = ', cminreal/(omegas(Nomega)/(2*Pi))
      write(*,*) '#intervalos por comprimento de onda (no modelo)= ',  ( cminreal/(omegas(Nomega)/(2*Pi)))/delta

      if(iprint.ne.0)then !se iprint =0 nao plotar o modelo de velocidades

	  call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloc',&
	    & 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
	    & wbox,Lz0/4.d0,Lx0/4.d0, 0.0d0,Lz0,0.0d0,Lx0,1)!modelo do pd sem PML

	  call systempsimage(rho,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'densi',&
	    & 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
	    & wbox,Lz0/4.d0,Lx0/4.d0, 0.0d0,Lz0,0.0d0,Lx0,1)!modelo do pd sem PML
      endif


      if(tsim.eq.1)then ! modelagem no domínio da frequencia

	  ! 	      write(*,*)'alocando P(Nx,Nz,Nomega,NsourceSnap)'
	  ! 	      allocate(P(Nx,Nz,Nomega,NsourceSnap))

	  write(*,*)'alocando P(Nx,Nz,NsourceSnap,gruposSR)'
	  allocate(P(Nx,Nz,NsourceSnap,gruposSR))
	  write(*,*)'alocando Ptreal(Nx,Nz,Ntsnaps,NsourceSnap,gruposSR)'
	  allocate(Ptreal(Nx,Nz,Ntsnaps,NsourceSnap,gruposSR))
	  ! 	      allocate(Ptimag(Nx,Nz,Ntsnaps,NsourceSnap))
	  ! 	      dtsnaps = (tf/dfloat(ntsnaps-1))

	  ! 	      ideltasnap = int(Nsource/(Nsourcesnap-1))


! 	  nttotal = 2*Nomega+1
	  dt=(tf/dfloat(nttotal))

	  do kw = 1,Nomega

	      write(*,*)'.... kw',kw
	      call PDMDFDF1(omegas(kw),0,0,0)!< iflagF=0 iflagG=0 iflagH=0

	      do igsr=1,gruposSR !loop nos grupos fonte receptor
		  ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))
		  do kss=1,NsourceSnap ! acumula resp no dom do tempo p cada freq resolvida para algumas fontes
		      ks=(kss-1)*ideltaSnap + 1
		      if(ks.gt.Nsource(igsr))ks=Nsource(igsr)
		      do it=1,ntsnaps
! 			    t = (it-1)*dtsnaps
			  t=(it-1)*nttotal/(ntsnaps-1)*dt ! para alguns ntsnaps instantes usados na transformada da fonte 
			  do j=1,Nz
			      do i=1,Nx
				  ! parte real do campo depois da tranformada inversa
				  Ptreal(i,j,it,kss,igsr) = Ptreal(i,j,it,kss,igsr) + real(          P(i,j,kss,igsr)/tf*exp(imag*omegas(kw)*t) + &  
										  & CONJG(P(i,j,kss,igsr))/tf*exp(imag*omegas(2*Nomega-(kw-1))*t))
				  ! parte imaginaria do campo (praticamente nula ja que f eh real)
! 				    Ptimag(i,j,it,kss) = Ptimag(i,j,it,kss) +aimag(  P(i,j,kss)*exp(imag*omegas(kw)*t) + &  
! 										    & CONJG(P(i,j,kss))*exp(imag*omegas(2*Nomega-(kw-1))*t)) 
			      enddo! em x
			  enddo! em z
		      enddo ! fim do loop nos instantes
		  enddo! nas fontes escolhidas para os snapshots
	      enddo ! loop os grupos fonte receptor


	  enddo! loop das frequencias


	  open(1,file='medidasfreq.txt')
	  write(1,*)'#Nreceptor Nfonte Nfrequencia Medida'
	  do kw=1,Nomega
	      do igsr=1,gruposSR
		  do ks=1,Nsource(igsr)
		      do kr=1,Nrec(igsr)
			  write(1,*)Preceiv(kr,ks,kw,igsr)
! 			  write(1,101)kr,ks,kw,Preceiv(kr,ks,kw)
		      enddo
		  enddo
	      enddo
	  enddo
	  close(1)

	  ! Faz a transformada inversa calculando a resp no dominio do tempo em cada receptor e gera sismograma
	  write(*,*)'FFT inversa e sismogramas'

! 	  nt = 2*Nomega+1
! 	  dt=(tf/dfloat(nt))


	  do igsr=1,gruposSR



	      allocate(PtrecR(nttotal,Nrec(igsr)))
! 	      allocate(PtrecI(nttotal,Nrec(igsr)))
! 	      ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))
! 
! 	      do kss=1,NsourceSnap
! 
! 		  ks=(kss-1)*ideltaSnap + 1
! 		  if(ks.gt.Nsource(igsr))ks=Nsource(igsr)

		  kss=1

		  do ks=1,Nsource(igsr)

		      do it=1,nttotal
			  t = (it-1)*dt
			  do kr=1,Nrec(igsr)
			      Ptrecaux = 0.0d0
			      do kw = 1,Nomega
				      Ptrecaux = Ptrecaux + Preceiv(kr,ks,kw,igsr)/tf*exp(imag*omegas(kw)*t)+&
					    & conjg(Preceiv(kr,ks,kw,igsr))/tf*exp(imag*omegas(2*Nomega-(kw-1))*t)
			      enddo
			      PtrecR(it,kr) = real(Ptrecaux)
! 			      PtrecI(it,kr) =aimag(Ptrecaux)

			  enddo
		      enddo

		      if(ks.lt.10)then
			  write(chks,'(i1)')ks
			  chks = '0000'//chks
		      elseif(ks.lt.100)then
			  write(chks,'(i2)')ks
			  chks = '000'//chks
		      elseif(ks.lt.1000)then
			  write(chks,'(i3)')ks
			  chks = '00'//chks
		      elseif(ks.lt.10000)then
			  write(chigsr,'(i4)')ks
			  chks = '0'//chks
		      elseif(ks.lt.100000)then
			  write(chks,'(i5)')ks
		      endif
		      if(igsr.lt.10)then
			  write(chigsr,'(i1)')igsr
			  chigsr = '0'//chigsr
		      elseif(igsr.lt.100)then
			  write(chigsr,'(i2)')igsr
		      endif

		      dreceiv=dfloat(ireceiv(2,igsr)-ireceiv(1,igsr))*delta ! distancia entre receptores
		      call systempswigp(PtrecR,nttotal,nrec,dt,dreceiv,&
				& 0.0d0,dfloat(ireceiv(1,igsr)-isource(ks,igsr))*delta,&
				& filesism//chigsr//chks,'      sismograma     ','      tempo (s)      ','     offset (m)      ',&
				& wbox,tf/5.0d0,dfloat(ireceiv(nrec(igsr),igsr)-ireceiv(1,igsr))*delta/5.d0,iprint)

		      if(ks.eq.kss)then
			  kss = kss + ideltasnap ! já atualiza a próxima fonte para a qual o sismograma será mantido 
			  if(kss.gt.Nsource(igsr))kss = Nsource(igsr)! corrige o número da fonte para plotagem do sismograma
		      else
			  call system('rm '//filesism//chigsr//chks//'.ps')! deletando alguns sismogramas para poupar espaço em disco
		      endif

! 		      if(ks.eq.30)then
! 			  kr=30
! 			  open(unit=10,file='traco-df.txt')
! 			  do it=1,nttotal
! 			      write(10,*)(it-1)*dt,PtrecR(it,kr)
! 			  enddo
! 			  close(10)
! 		      endif

		  enddo ! loop nas fontes

! 	      enddo! loop das fontes nsourcesnap
	      

	      deallocate(PtrecR)
! 	      deallocate(PtrecI)

	  enddo ! fim do loop dos grupos fonte receptor
	  !  fim da geração dos sismogramas


	  do igsr=1,gruposSR
	  ! a partir da resp ja no dom do tempo  gera os snapshots
	      write(*,*)'snapshots'
	      ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))
	      if(igsr.lt.10)then
		  write(chigsr,'(i1)')igsr
		  chigsr = '0'//chigsr
	      elseif(igsr.lt.100)then
		  write(chigsr,'(i2)')igsr
	      endif


	      do kss=1,NsourceSnap

		  ks=(kss-1)*ideltaSnap + 1
		  if(ks.gt.Nsource(igsr))ks=Nsource(igsr)

		  if(ks.lt.10)then
		      write(chks,'(i1)')ks
		      chks = '0000'//chks
		  elseif(ks.lt.100)then
		      write(chks,'(i2)')ks
		      chks = '000'//chks
		  elseif(ks.lt.1000)then
		      write(chks,'(i3)')ks
		      chks = '00'//chks
		  elseif(ks.lt.10000)then
		      write(chks,'(i4)')ks
		      chks = '0'//chks
		  elseif(ks.lt.100000)then
		      write(chks,'(i5)')ks
		  endif


		  filename4='snaps'//chks

		  open(1,file=filename4//'.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx*Nz)
		  do it=1,ntsnaps
		      write(1,rec=it)((Ptreal(i,j,it,kss,igsr),j=1,Nz),i=1,Nx)
		  enddo
		  close(1)      

		  dtsnaps = tf/dfloat(ntsnaps-1)


		  call systempsmovie(nx,nz,delta,0.00d0,-dfloat(Npml-1)*delta,0.0d0,dtsnaps,&
					      &'        campo       ','  profundidade (m)  ',&
					      &'   distancia (m)    ',wbox,Lz0/4.d0, Lx0/4.d0 ,  &
					      & ideltasnap, ntsnaps,filename4  )
		  !   fim da geracao dos snapshots para algumas fontes

	      enddo! loop das fontes
	  enddo! loop dos grupos de carregamento

	  deallocate(P)
	  deallocate(Ptreal)
	  !deallocate(Ptimag)


      elseif(tsim.eq.2)then ! modelagem no domínio do tempo

	  write(*,*)'nttotal', nttotal
	  nttotalIN = nttotal
	  ntpsIN = ntps
	  deltatIN = 1.0d0/dfloat(ntps)


	  deltat = 1.0d0/dfloat(ntps)

	  ! verificacao da discretizacao
	  cminmodel = c(1,1)
	  cmaxmodel = c(1,1)
	  do i=1,nx
	      do j=1,nz
		  if(c(i,j).lt.cminmodel)then
		      cminmodel=c(i,j)
		  endif
		  if(c(i,j).gt.cmaxmodel)then
		      cmaxmodel=c(i,j)
		  endif
	      enddo
	  enddo
	  deltamax=cminmodel/(5.*fcorte)
	  deltatmax=deltamax/(5.*cmaxmodel)

	  write(*,*) 'simulacao no domínio do tempo'
	  write(*,*) 'delta T maximo =',deltatmax
	  write(*,*) 'delta maximo =',deltamax
	  write(*,*) 'delta T =',deltat
	  write(*,*) 'delta =',delta
	  if(delta.gt.deltamax)then
	      write(*,*) 'delta usado maior que valor maximo permitido'
	  endif
	  if(deltat.gt.deltatmax)then
	      write(*,*) 'deltaT usado maior que valor maximo permitido... corrigindo...'
	      ntps = int(1.0d0/deltaTmax)+1
	      deltaT = 1.0d0/dfloat(ntps)
	      nttotal = int(tf*ntps)
	      write(*,*) 'delta T =',deltat
	      write(*,*) 'ntps =',ntps
	      write(*,*) 'nttotal =',nttotal
	  endif

	  allocate(St(nttotal),stat=ierr)
	  call rickerwavelet(tf,nttotal,fcorte,fpico,ampS,tp,St)


	  !alocando as matrizes auxiliares
	  allocate(Caux(nx,nz),stat=ierr)
	  if(ierr.eq.0)then
	      write(*,*)'.........Ok'
	  else
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif
	  allocate(Aaux(nx,nz),stat=ierr)
	  if(ierr.eq.0)then
	      write(*,*)'.........Ok'
	  else
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif

	  !inicializando as matrizes auxiliares
	  do j=1,nz
	      do i=1,nx
		  Caux(i,j) =-(c(i,j) *deltat / delta)**2 / 12.d0
		  Aaux(i,j) = c(i,j) * deltat / delta
	      enddo
	  enddo

	  ! alocando as matrizes solução
	  allocate(P1(nx,nz),stat=ierr)
	  if(ierr.eq.0)then
	      write(*,*)'.........Ok'
	  else
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif
	  allocate(P2(nx,nz),stat=ierr)
	  if(ierr.eq.0)then
	      write(*,*)'.........Ok'
	  else
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif
	  allocate(P3(nx,nz),stat=ierr)
	  if(ierr.eq.0)then
	      write(*,*)'.........Ok'
	  else
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif



	  do igsr=1,gruposSR ! Loop no número de grupos fonte/receptor

	      ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))! de qtas em quantas fontes será feito o snapshot para este grupo fontes/receptores
	      dtsnaps = tf/dfloat(ntsnaps-1)
	      if(igsr.lt.10)then
		  write(chigsr,'(i1)')igsr
		  chigsr = '0'//chigsr
	      elseif(igsr.lt.100)then
		  write(chigsr,'(i2)')igsr
	      endif
	      kss = 1
       
	      do ks =1,Nsource(igsr) ! Loop nas fontes

		  if(ks.lt.10)then !salvar ks no formato para sismogramas e snapshots
		      write(chks,'(i1)')ks
		      chks = '0000'//chks
		  elseif(ks.lt.100)then
		      write(chks,'(i2)')ks
		      chks = '000'//chks
		  elseif(ks.lt.1000)then
		      write(chks,'(i3)')ks
		      chks = '00'//chks
		  elseif(ks.lt.10000)then
		      write(chks,'(i4)')ks
		      chks = '0'//chks
		  elseif(ks.lt.100000)then
		      write(chks,'(i5)')ks
		  endif


		  if(ks.eq.kss)then
		      plotsnap=1 !<vai plotar o snapshot

		      filename4='snaps'//chks ! abre o arquivo para plotagem do campo
		      open(1,file=filename4//'.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx*Nz)

		      kss = kss + ideltasnap ! já atualiza a próxima fonte para a qual será salvo o snapshot
		      if(kss.gt.Nsource(igsr)) kss = Nsource(igsr)! corrige o número da fonte para o próximo snapshot caso este número seja maior que o número de fontes
		  else
		      plotsnap=0 !<nao vai plotar o snapshot
		  endif


		  allocate(PtrecR(nttotal,Nrec(igsr)))



		  ! para a fonte atual, calcula o campo em todo o domínio
		  call PDMDFDT1(nttotal,ntps,St,nx,nz,delta,deltat,c,rho,&
		      &  isource(ks,igsr),jsource(ks,igsr),Nrec(igsr),ireceiv(1,igsr),jreceiv(1,igsr),&
		      & P1,P2,P3,Aaux,Caux,plotsnap,ntsnaps,PtrecR)

		  if(plotsnap.eq.1)then!plotando snapshots, se for o caso
		      close(1) !fecha o arquivo bin que foi escrito dentro da modelagem
		      call systempsmovie(nx,nz,delta,0.00d0,-dfloat(Npml-1)*delta,0.0d0,dtsnaps,&
					      &'        campo       ','  profundidade (m)  ',&
					      &'   distancia (m)    ',wbox,Lz0/4.d0, Lx0/4.d0 ,  &
					      & ideltasnap, ntsnaps,filename4  )
		      !   fim da geracao dos snapshots para algumas fontes
		  endif







		  !plotando sismograma
		  if(nttotal.gt.nttotalIN)then !verifica se os sismogramas foram ferados com mais instantes de tempo que a taxa de amostragem desejada

		      allocate(PtrecRaux(nttotal,Nrec(igsr)))!fazer uma copia dos sismogramas com mais instantes de tempo
		      do kr=1,Nrec(igsr)
			  do it=1,nttotal
			      PtrecRaux(it,kr) = PtrecR(it,kr)
			  enddo
		      enddo

		      deallocate(PtrecR)
		      allocate(PtrecR(nttotalIN,Nrec(igsr)))! alocando os sismogramas em uma matriz menor


		      do itIN=1,nttotalIN
			  tIN = (itIN-1)*1.0d0/dfloat(ntpsIN)

			  do it=1,nttotal-1
			      t1=(it  -1)*1.0d0/dfloat(ntps)
			      t2=(it+1-1)*1.0d0/dfloat(ntps)
			      it1 = it
			      it2 = it+1
			      if(tIN.ge.t1 .and. tIN.le.t2)exit
			  enddo

			  do kr = 1, Nrec(igsr)
			      PtrecR(itIN,kr) = PtrecRaux(it1,kr)+ (PtrecRaux(it2,kr)-PtrecRaux(it1,kr))*(tIN-t1)/(t2-t1)
			  enddo

		      enddo


		      deallocate(PtrecRaux)

		  endif




		  write(*,*) 'plotando sismograma ',filesism,chigsr,chks,'.bin'
		  dreceiv=dfloat(ireceiv(2,igsr)-ireceiv(1,igsr))*delta ! distancia entre receptores
		  call systempswigp(PtrecR,nttotalIN,nrec(igsr),deltatIN,dreceiv,&
			& 0.0d0,dfloat(ireceiv(1,igsr)-isource(ks,igsr))*delta,&
			& filesism//chigsr//chks,'      sismograma     ','      tempo (s)      ','     offset (m)      ',&
			& wbox,tf/5.0d0,dfloat(ireceiv(nrec(igsr),igsr)-ireceiv(1,igsr))*delta/5.d0,iprint)

		  if(plotsnap.eq.0)then ! se não vai plotar snapshot, aproveita para deletar s sismogramas
		      call system('rm '//filesism//chigsr//chks//'.ps')! deletando alguns sismogramas para poupar espaço em disco
		  endif

! 		  if(ks.eq.30)then
! 		      kr=30
! 		      open(unit=10,file='traco-dt.txt')
! 		      do it=1,nttotal
! 			  write(10,*)dfloat(it-1)*tf/dfloat(nttotal),PtrecR(it,kr)
! 		      enddo
! 		      close(10)
! 		  endif





		  deallocate(PtrecR)


	      enddo ! fechamento do Loop nas fontes
	  enddo ! fechamento do  Loop no número de grupos fonte/receptor
   

	  deallocate(P1)
	  deallocate(P2)
	  deallocate(P3)

	  deallocate(Caux)
	  deallocate(Aaux)	

	  nttotal = nttotalIN!caso tenha sido feita a simulação no tempo com passo de tempo diferente do determinado pela entrada, retorna aos valores da entrada
	  ntps = ntpsIN
	  deltat = deltatIN      

      endif! fim do loop de escolha do tipo de simulação

      call system('tar --remove-files -czf '//filesism//'.tgz '//filesism//'*.bin')! adiciona ao arquivo compactado deletando os arquivos



      return

      end subroutine simulacao